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Abstract 

Phase fitting has been extensively used during the last years to improve the 
behaviour of numerical integrators on oscillatory problems. In this work, the 
benefits of the phase fitting technique are embedded in discrete Lagrangian 
integrators. The results show improved accuracy and total energy behaviour 
in Hamiltonian systems. Numerical tests on the long term integration (10 5 
periods) of the 2-body problem with eccentricity even up to 0.95 show the 
efficiency of the proposed approach. Finally, based on a geometrical eval- 
uation of the frequency of the problem, a new technique for adaptive error 
control is presented. 
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1. Introduction 



In the field of numerical integration, methods specially tuned on os- 
cillating functions, are of great practical importance. Such methods are 
needed in various branches of natural sciences, particularly in physics, since 
a lot of physical phenomena exh ibit a prono u nced oscillatory behaviour . 
For a review of such methods see ISimosI (120001); |L. Ixaru and Daek- (1997); 



G. Vanden Berghe and Hecke 
L. Or. Ixaru and Meverl (l200. c 



(19 99h: IG. Vanden Berghe and Daeld ((20011); 



Daele and Berghd (120071 ) and references there 



in as well as the book lL.Gr. Ixaru! (120041 ) . 
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For problems having highly oscillatory solutions, standard methods with 
unspecialised use can require a huge number of steps to track the oscilla- 
tions. One way to obtain a more efficient integration process is to construct 
numerical methods with an increased algebraic order, although the simple 
implementation of high algebraic order methods may cause severa l problems 
(for example, the existence of parasitic solutions iQuinlanl (119991 ) ) . On the 
other hand, there are some special techniques for optimising numerical meth- 
ods. Trigonometrical fitting and phase-fitting are some of them, producing 
methods with variable coefficients, which depend on v — uh, where u is the 
dominant frequency of the problem and h is the step length of integration. 
This technique is kno wn as exponentia l (or trigonom etric if \i — ioj) fitting 
and has a long history iGautschil (119611 ). iLychd (Il972l ). An important prop- 
erty of exponential fitted algorithms is that they tend to the classical ones 
when the involved frequencies tend to zero, a fact which allows to say that 
exponential fitting represents a natural extension of the classical polynomial 
fitting. The examination of the conve rgence of exp onential fitted multistep 
methods is included i n Lyches theory iLvchd (119721 ). The general theory is 
presented in detail in lL.Gr. Ixarul (120041 ) . Furthermore, considering the ac- 
curacy of a method when solving oscillatory problems, it is more appropriate 
to work with the phase-lag, rather than its usual primary local truncatio n 
error. We mention the pioneering paper of Brusa and Nigro lL. Brusal (119801 ). 
in which the phase-lag property was introduced. This is actually another 
type of a truncation error, i.e. the angle between the analytical solution and 
the numerical solution. A significant application of the phase or exponential 
fitting is on the construction of sympl ectic methods for oscillatory problems 
encounterd in physics and chemistry ( T. Monovasilis ( 2005 . 2006I )). 

Another approach to oscillatory and especially Hamiltonian systems is 
the theory of disc r ete va r iational mecha n ics, wh i ch wa s set up in the 1960s 
Jordan and Polak ( 1964 ); Cadzow ( 197Cll ); Logan ( 1973 ) and then it was pro- 
posed in the optimal control literature. It then motivated a lot of authors 
and soon the discrete Euler-Lagrange equations were formulated and the 
first integrators in the discrete calculus of variation and further the multi- 
freedom and higher-order problems were studied. Afterwards, the canonical 
structure and symmetries for discrete system s were obtained, and Noether's 
theorem to the discrete case was extended iMaedal (119801. Il98 lh . Finally, 



the time as a discrete dynamical variable was regarded iLed (119831 ). A de- 
tailed de scription of the essential p ro perties of variational integr a tors c a n be 
found in iMarsden and Westl (l200lh : fj.E. Marsden and Shkollerl (Il998h : [Lee 



(119831 ) . One of the most important properties of variational integrators is 
that since the discrete Lagrangian is an approximation of a continuous La- 
grangian function, the obtained numerical integrator inherits some of the 
geometric properties of the continuous Lagrangian (such as symplecticity, 
momentum preservation). 

In the present work, the benefits of the two approach are combined in or- 
der to construct discrete Lagrangian integrators with phase fitting. To obtain 
this, we have adopted a test Lagrangian problem (similar to the test ODE 
in the phase fitting) which is the harmonic oscillator with given frequency uj. 
Then, we construct discrete variational schemes that solve exactly the test 
Lagrangian. The application of the method to a general Lagrangian needs 
the determination of the frequency u at every step of the integration. The 
method is applied to the 2-body problem with eccentricity up to 0.95. The 
results exhibit an improved behaviour of the calculated solutions, especially 
in the case of total energy of the integrated systems. 



2. Discrete Variational Mechanics 

The well known least action principle of the continuous Lagrange - Hamil- 
ton Dynamics can be used as a guiding principle to derive discrete integra- 
tors. Following the steps of the derivation of Euler-Lagrange equations in 
the continuous time Lagrangian dynamics, one can derive the discrete time 
Euler-Lagrange equations. For this purpose, one considers positions q and 
qi and a time step h G R, in order to replace the parameters of position q 
and velocity q in the continuous time Lagrangian L(q,q,t). Then, by con- 
sidering the variable h as a very small (positive) number, the positions qo 
and qi could be thought of as being two points on a curve (trajectory of the 
mechanical system) at time h apart. Under these assumptions, the following 
approximations hold: 

Qo « 9(0) , qi w q(h) , 

and a function L d (q , q X) h) could be defined known as a discrete Lagrangian 
function. 

Many authors assume such functions to approximate the action integral 
along the curve segment between q and qi, i.e. 

Ld(qo,qi,h)= L(q,q,t)dt (1) 
Jo 
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Furthermore, one may consider the very simple approxi mation for this inte 



Krai g iven on the basis of the rectangle rule described in iMarsden and West 
(120011 ). According to this rule, the integral f Q Ldt could be approximated 



by the product of the time-interval h times the value of the integrand L ob- 
tained with the velocity q replaced by the approximation (qi — qo)/h: The 
next step is to consider a discrete curve defined by the set of points {q k } k=0 , 
and calculate the discrete action along this sequence by summing the discrete 
Lagrangian of the form Ld(q k , qk+i, h) defined for each adjacent pair of points 

(qk, qk+i)- 

Following the case of the continuous dynamics, we compute variations 
of this action sum with the boundary points go an d qN held fixed. Briefly, 
discretization of the action functional leads to the concept of an action sum 



ra-1 



Sdipld) = ^2 L d(lk-h Qk), Id = Qn-l) G Q n (2) 



k=l 



where La : Q x Q — > R is an approximation of L called the discrete La- 
grangian. Hence, in the discrete setting the correspondence to the velocity 
phase space TQ is Q x Q. An intuitive motivation for this is that two points 
close to each other correspond approximately to the same information as one 
point and a velocity vector. The discrete Hamilton's principle states that if 
7d is a motion of the discrete mechanical system then it extremizes the ac- 
tion sum, i. e., 5Sd = 0. By differentiation and rearranging of the terms and 
having in mind that both q and q^ are fixed, the discrete Euler-Lagrange 
(DEL) equation is obtained: 

D 2 L d (qk-i, qk, h) + D x L d {qk, qk+i, h) = (3) 

where the notation DiL d indicates the slot derivative with respect to the 
argument of Z^. 

In a position-momentum form the discrete Euler-Lagrange equations ([3]) 
can be defined by the equations below 

p k = -D 1 L d (q k ,q k+1 ,h) 
p k+1 = D 2 L d (q k , q k+1 , h) (4) 

3. Phase-fitted Discrete Lagrangian Integrators 

Summarising the phase fitting technique, we consider for simplicity only 
first order differential equations, although the same results can be easily 
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obtained for second order equations too. Consider the test problem 



^- = iu y(t), y(0) = l (5) 

with exact solution 

y(t) = e*"* (6) 

where is a non-negative real value. Let &(h) be a numerical map which 
when it is applied to a set of known past values, it produces a numerical 
estimation of y(t + h). If we assume that all past values are known exactly, 
then the numerical estimation y(t + h) of y(t + h) will be 

y(t + h) = a{u h) ■ e ^ot+H^oh)) ( 7 ) 

while the exact solution is e l ^°' +aJ0?t ). Then the ratio of the estimated to the 
exact solution is 

= y(t±h) = a ^ h ) e -i^oh-^ h)) (8) 

y(t + h) v ; K J 

In the above equation ([8]), the term a(u>oh) is called the amplification error, 
while the term 1{ujqK) = uj h — <f)(u)oh) is called the phase lag of the numerical 
map. In the case that a(uioh) = 1 and 1{ojqK) = 0, we say that the numerical 
map is exponentially fitted at the frequency ujq and at the step size h. 
The technique of phase fitting can now be considered as the vanishing or 
minimisation of the phase lag. 

Consider now the discrete Lagrangian Ld(qk,qk+i,h) (qk corresponds to 
time tk and qu+i to time t^+i = tfc + h) and a set of s intermediate points 
g J with qi = q (t^ + c'h). The role of the number of intermidiate points will 
be discussed later. Assuming that c 1 = and c s = 1 we always have q 1 = q^ 
and q s = qk+i- Then we can approximate L d with the quadrature 

s 

L d (q k , q k+1 , h) = h^Wj ■ L (q(t k + c j h),q(t k + c j h),c j h) (9) 

For maximal algebraic order it is easily proved that the following conditions 
must hold: 

E^M' = rn;> i = o,i,.. (io) 
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Then, we can approximate intermediate points and their derivatives with 

qi = b>q k + b>q k+ i , , 

Consider now the test Lagrangian (harmonic oscillator) similar to the test 
equation (JSJ) 

u = \e - (12) 

Then, applying the above assumptions in Eq. ([3]) we get 



e; =1 ^((^) 2 +(^) 2 -« 2 ((^) 2 +(^) 



where u = uh. Since the exact solution of Eq. (fl2|) is 

= Ae ibA + Be~ iLJt (14) 
and we want to force our method to solve exactly Eq. we get 
V = cos {c?u) — cosu s i n {r?u\ 

V 



cosu 
sinu 

sin (c?u) 



sinu 

cosu 



B 3 = —usin (c?u) — u— — cos (c^u) 



sinu 



cos (c J u) 



B 3 = u ^ '- (15) 

sinu 

4. Frequency Evaluation and Error Control 

The final step to our method is to evaluate the frequency of the problem. 
For this purpose we use the curvature of the solution and the concept of 
the oscullating circle. Consider a planar curve C and a point P on this 
curve. If r(t) is a parametrised representation of the curve, then we define 
the curvature at a point P as 

, , . r(t) x r(t) , . 

k{t) = IMF (16) 
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Then, there is a circle with radius 

R 



\k(t)\ 



(17) 



which locally approximates the curve at point r(t) and is called the oscul- 
lating circle. Since, the velocity of a point running on top of the curve C is 
|r(t)|, at a small time step h the point will rotate an angle equal to 

|r(t)| , |r x f | 
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R 



h 



leading us to a frequency selection 



r x r 
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(18) 



(19) 



Consider now the regular parametrisation of the curve r (the one that uses 
the curve length as the free parameter). Let s the curve length. Then the 
curve u of the centres of the oscullating circles is given by 

1 



us 



r is) 



k(s 



(20) 



Then, 



(21) 



where H(s) is the first normal vector of the curve r at point r(s) 

du(s) dk(s)/ds^^ , 
ds k(s) 2 

which means that the centre of the oscullating circle is moving to a direction 
normal to the curve r with velocity 

dk(s)/ds 



k(s) 



(22) 



Assuming now a small displacement on curve r, it can be easily proved that 
the distance that is covered by the centre of the oscullating circle is smaller 
than the difference of their radius which means that the one circle is entirely 
inside or outside of the other (depending on the variation of the curvature). 
Thus, the error in position is bounded by 

1 1 

h k 2 

where k±, k 2 are the curvatures of the curve r at the adjacent points and h 
is the time step. Figure §I§ depicts this result. Using Eq. (123]) we can adap- 
tively control the time step of the integration, keeping the local truncation 
error within desired bounds. 



- v h 



(23) 
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Figure 1: The oscullating circles at two neighbour points P\ and P 2 are shown. The first 
circle has its centre at point O and radius Rl and the second its centre at point O' with 
radius R2. If P 2 is the estimated point and the time step h is small enough, then the 
distance between P 2 and P 2 is of the order of the absolute difference of the two radii 
minus the distance between O and O' . 

5. Numerical Test for the 2-body Problem 

We now turn to the study of two objects interacting through a central 
force. The most famous example of this type, is the Kepler problem (also 
called the two-body problem) that describes the motion of two bodies which 
attract each other. In the solar system the gravitational interaction between 
two bodies leads to the elliptic orbits of planets and the hyperbolic orbits of 
comets. 

If we choose one of the bodies as the centre of our coordinate system, the 
motion will stay in a plane. Denoting the position of the second body by 
q = (qi, q2) T , the Lagrangian of the system takes the form (assuming masses 
and gravitational constant equal to 1) 




(24) 



The initial conditions are taken 




(25) 
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where e is the eccentricity of the orbit. In order to check the efficiency of 
the proposed algorithm, we shall consider only high eccentricity (e = 0.95). 
Figure ([2]) compar e s the proposed method with the method described in 



Marsden and West ( 200ll ) (the two methods have the same algebraic order). 



The results here are obtained as follows: First, the phase fitted method is 
applied for one period and for a given tolerance in positio n calculation. Then, 



the en ergy tolerance is calculated and the method of iMarsden and West 



(120011 ) is applied using a variable step length in order to obtain the same 
energy tolerance. It is clear that the phase fitting decreases close to one 
third the number of integration steps to obtain the same accuracy. 

Finally, in figure [3] the total energy is plotted as a function of time for 10 5 
periods as well as the error in position (the distance between the calculated 
and the exact points). It is clear that the method keeps both energy and 
position error in stable limits although there is a relative increased error in 
energy at the perihelion. This can be explained by considering the calculated 
from equation Eq. ffl9l) frequency of the problem(see Fig. H|). The frequency 
is smooth enough almost everywhere except at the perihelion where it changes 
rapidly. This means that the assumption that the frequency is constant 
during an integration step, applies everywhere except at the perihelion and 
this is the reason for the observed increase in the total energy error. This 
undesirable effect can be handled by increasing the algebraic order of the 
method. 



6. Increasing the algebraic order 

Until now, he have not mention anything about the role of the number 
of intermediate points used to calculate the discrete Lagrangian. This can 
be used to increase the algebraic order of the method. Consider Eq. (fTTj) 
modified as 

q 3 = Vq k + Vq k+ i + bq 3 , , 

q 3 = \ (B 3 q k + B 3 q k+l ) + 5q 3 { ° J 

where the corrections Sq 3 are free parameters and the corresponding correc- 
tions Sq 3 for the derivatives cane be calculated 

k=l 
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Absolut Error in Total Energy 



Figure 2: The number of integration steps needed to obtain the sa me accuracy in total 
energy for the proposed method (•) and for the method described in Marsden and West] 
lj200l( i (0). 




Figure 3: The total energy (solid line) end the error in position (dashed line) as a function 
of time for eccentricity 0.95 and for 10 5 periods. 
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6.28 12.56 



25.12 31.4 37.6 

Time 



43.96 50.24 56.52 



Figure 4: The calculated frequency of the problem is shown for the 10 first periods. 

where the coefficients a 3 k can be easily calculated for maximal algebraic order 

as 







i=i 

±4W = »(<*) 



n-l 



n= 1,2, 



(28) 



Now, the discrete Lagrangian, beyond (qk,qk+i) depends also on 8q> . The 
system of Eq. (J4j) now is enriched with the equations (since we want the 
discrete Lagrangian to be stationary) 

dL d (q k , q k+1 , h) 



d5qi 



0, j = l,2,. ..,8 



(29) 



This t echnique is similar to those described in lLeold (120051 ) and lL. Kharevych and Desbrun 
fl2006h . 



7. Conclusions 

It has been shown in this work, that the technique of phase fitting, when 
it is embedded in discrete Lagrangian integrators, improves the accuracy and 
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the energy behaviour of the numerical method. Following the classical ap- 
plication of the phase fitting technique, the discrete Lagrangian integrator 
is forced to solve exactly the test Lagrangian of harmonic oscillator with a 
given self frequency. The coefficients of the resulting integrator, depend on 
the frequency of the problem at each integration step. Geometrical consid- 
eration lead us to a new frequency evaluation depending on the curvature 
of the solution and on the principle of the oscullating circle. Furthermore, 
the resulting analysis gave us a new error control technique. Application of 
the new method to the well known 2-body problem decreased the number 
of integration steps needed to obtain the same accuracy with the classical 
discrete Lagrangian of the same algebraic order to less than one third in high 
eccentricity equal to 0.95. Moreover, the new method exhibits improved en- 
ergy behaviour for long term integration (10 5 periods of the 2-body problem 
with eccentricity equal to 0.95). Finally, a simple method is proposed to 
increase the algebraic order of the new method. 

8. Acknowledgement 

This paper is part of the 03ED51 research project, implemented within 
the framework of the " Reinforcement Programme of Human Research Man- 
power" (PENED) and co-financed by National and Community Funds (25% 
from the Greek Ministry of Development-General Secretariat of Research and 
Technology and 75% from E.U. -European Social Fund). 

References 

Cadzow, J. A., 1970. Discrete calculus of variation. Internat. J. Control 11, 
393-407. 

Daele, M. V., Berghe, G. V., 2007. Geometric numerical integration by means 
of exponentially fitted methods. APNUM 57, 415-435. 

G. Vanden Berghe, L. G. L, Daele, M. V., 2001. Optimal implicit 
exponentially-fitted runge-kutta methods. Computer Physics Comm 150, 
346-357. 

G. Vanden Berghe, H. De Meyer, M. V. D., Hecke, T. V., 1999. 
Exponentially-fitted explicit runge-kutta methods. Computer Physics 
Comm 123, 7-15. 



12 



Gautschi, W., 1961. Numerical integration of ordinary differential equations 
based on trigonometric polynomials. Numer. Math. 3, 381397. 

J.E. Marsden, G. P., Shkoller, S., 1998. Multisymplectic geometry, variational 
integrators and non-linear pdes. Comm. Math. Phys. 199, 351-395. 

Jordan, W., Polak, E., 1964. Theory of a class of discrete optimal control 
systems. J.Eletron. Control 17, 697-711. 

L. Brusa, L. N., 1980. A one-step method for direct integration of structural 
dynamic equations. Int. J. Num. Methods Engrg. 15, 685699. 

L. Gr. Ixaru, G. V. B., Meyer, H. D., 2003. Exponentially fitted variable 
two-step bdf algorithm for first order odes. Computer Physics Comm 150, 
116-128. 

L. Ixaru, G. Vanden Berghe, H. D. M., Daele, M. V., 1997. Four-step 
exponential-fitted methods for nonlinear physical problems. Computer 
Physics Comm 100, 56-70. 

L. Kharevych, W.Y. Tong, E. K. J. M. P. S., Desbrun, M., 2006. Geometric, 
variational integrators for computer animation. Eutographics/ACM SIG- 
GRAPH Symposium on Computer Animation. 

Lee, T. D., 1983. Can time be a discrete dynamical variable? Phys. Lett. B 
122, 217-220. 

Leok, M., 2005. Generalized galerkin variational integrators. 
|arxiv:math/ 0508360 . 

L.Gr. Ixaru, G. V. B., 2004. Kluwer Academic Publishers, Dor- 
drecht / Boston / London. 

Logan, J. A., 1973. First integrators in the discrete calculus of variation. 
Aequationes Mathematicae 9, 210-220. 

Lyche, T., 1972. Chebyshevian multistep methods for ordinary differential 
equations. Num. Math. 19, 65-75. 

Maeda, S., 1980. Canonical structure and symmetries for discrete systems. 
Math. Japonica 25, 405-420. 



13 



Maeda, S., 1981. Extension of discrete noether's theorem. Math. Japonica 
26, 85-90. 

Marsden, J., West, M., 2001. Discrete mechanics and variational integrators. 
Acta Num. 10, 357-514. 

Quinlan, G., 1999. Resonances and instabilities in symmetric multistep meth- 
ods, preprint arXiv [astro- ph / 990 1 136 . 

Simos, T., 2000. Specialist Periodical Reports. The Royal Society of Chem- 
istry, Cambridge. 

T. Monovasilis, Z. Kalogiratou, T. S., 2005. Exponentially fitted symplec- 
tic methods for the numerical integration of the Schrodinger equation. J. 
Math. Chem 37 (3), 263-270. 

T. Monovasilis, Z. Kalogiratou, T. S., 2006. Trigonometrically fitted and 
exponentially fitted symplectic methods for the numerical integration of 
the Schrodinger equation. J. Math. Chem 40 (3), 257-267. 



14 



